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Abstract 

An efficient method, using equivalent plate model, is developed for studying the static and vibration 
analyses of general built-up wing structures composed of skins, spars, and ribs. The model includes the 
transverse shear effects by treating the built-up wing as a plate following the Reissner-Mindlin theory, the 
so-called First-order Shear Deformation Theory (FSDT). The Ritz method is used with the Legendre 
polynomials being employed as the trial functions. This is in contrast to previous equivalent plate model 
methods which have used simple polynomials, known to be prone to numerical ill-conditioning, as the trial 
functions. The present developments are evaluated by comparing the results with those obtained using 
MSC/NASTRAN, for a set of examples. These examples are: (i) free-vibration analysis of a clamped 
trapezoidal plate with (a) uniform thickness, and (b) non-uniform thickness varying as an airfoil, (ii)free- 
vibration and static analyses (including skin stress distribution) of a general built-up wing, and (iii)free- 
vibration and static analyses of a swept-back box wing. The results obtained by the present equivalent plate 
model are in good agreement with those obtained by the finite element method. 
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Introduction 


For structural analysis, the Finite Element Analysis (FEA) is widely used because of its generality, 
versatility and reliability. FEA is also the method of choice in situations where detailed results in the 
vicinity of local discontinuities (holes, abrupt dimension variations etc.) are needed. This is accomplished 
by refinement of the mesh near the zone of interest. But a wide application of detailed FEA at the late 
conceptual design stage or in the early preliminary design stage still faces some major obstacles. First, the 
preparation time for a FEA model data may be prohibitive, especially when there is little carry-over from 
design to design. Second, for complex structures, a detailed FEA needs huge amount of CPU time and 
computation capacity which makes the cost soar. 

In view of this situation, often equivalent continuum models are used to simulate complex structures 
for the purpose of obtaining global solutions in the early design stages. This idea is reasonable as long as 
the complex structure behaves physically in a close manner to the continuum model used and only global 
quantities of the response are of concern. For example, significant work exists on using beam or plate 
models to simulate repetitive lattice structures 1-7 . 

In the area of analyzing aerospace wing structures, a number of studies have been conducted on using 
equivalent beam models to represent simple box-wings composed of laminated or anisotropic 
materials 8 " 10 , and they have yielded accurate results for the specific problems studied. But from a realistic 
point of view, using an equivalent plate to model a wing structure is more promising, since a wing, 
especially the one that has a low aspect ratio, is likely to behave more as a plate than as a beam. 

There does exist a considerable body of work on the static or dynamic behaviors of all kinds of plates. 
A thorough description of literature on the study of plates was given by Lovejoy and Kapania 1 1 , where 
more than 300 references has been listed about all kinds of thin, thick, laminated or composite plates of 
trapezoidal shapes. One way of classifying existing methods for the solution of plates is according to the 
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deformation theory used, namely: the Classical Plate Theory (CPT), the First-order Shear Deformation 
Theory (FSDT), or the Higher-order Shear Deformation Theory (HSDT) etc. The CPT is based on the 
Kirchhoff-Love hypothesis, that is, a straight line normal to the plate middle surface remains straight and 
normal during the deformation process. This group of theories work well for truly thin isotropic plates, but 
for thick isotropic plates and for thin laminated plates they tend to overestimate the stiffness of the plate 
since the effects of through-the-thickness shear deformation are ignored 12 ' 13 . The FSDT is based on the 
Reissner-Mindlin model 14,15 , where the constraint that a normal to the mid-surface remains normal to the 
mid-surface after deformation is relaxed, and a uniform transverse shear strain is allowed. The FSDT is the 
most widely used theory for thick and anisotropic laminated plates owing to its simplicity and its low 
requirement for computation capacity. For more accurate results or more realistic local distributions of the 
transverse strain and stress, one should use the HSDT 16 , or the CFSDT (Consistent First-order Shear 
Deformation Theory) proposed by Knight and Qi 17 . 

Methods of solving the CPT, FSDT or HSDT mainly include finite element, Galerkin, and Rayleigh- 
Ritz methods 11 . In the context of using equivalent plate to represent the behaviors of wing structures at the 
conceptual stage at least, it is obvious that while the computationally costly finite element method is to be 
avoided, the Rayleigh-Ritz method is attractive. 

There have been several studies using equivalent plate models to model wing structures. Giles 1819 
developed a Ritz method based approach, which considers an aircraft wing as being formed by a series of 
equivalent trapezoidal segments, and represents the true internal structure of aircraft wings in the 
polynomial power form. In Giles 18 the CPT was used, but this shortcoming was removed subsequently 19 . 
Tizzi 20 presented a method whose many aspects are similar to that of Giles. In Tizzi’s work several 
trapezoidal segments in different planes can be considered, but the internal parts of wing structures (spars, 
ribs, etc.) were not considered. Livne 21 formulated the FSDT to be used for modeling solid plates as well 
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as typical wing box structures made of cover skins and an array of spars and ribs based on simple- 
polynomial trial functions, which are known to be prone to numerical ill-conditioning problems. Livne and 
Navarro then further developed the method to deal with nonlinear problems of wing box structures 22 . 

The present study is an extension of the previous works of Kapania and Singhvi 23,24 , Kapania and 
Lovejoy 25,26,1 1 and Cortial 27 , who all used the Rayleigh-Ritz method with the Chebyshev polynomials as 
the trial functions, and applied the Lagrange’s equations to obtain the stiffness and mass matrices. In 
Kapania and Singhvi 23 24 , the CPT was used to solve generally laminated trapezoidal plates, while in 
Kapania and Lovejoy 25,26,1 1 , the FSDT was used. In all these studies, only uniform plates were considered. 
In Cortial 27 , efforts were made to use the method of Kapania and Lovejoy 25,26,11 to calculate natural 
frequencies of box-wing structures, but an assumption of constant wing thickness makes it difficult to 
apply the method to general wing structures. 

In this paper, it is assumed that the wing plan-form is trapezoidal, and the wing structure is composed 
of skins, spars and ribs. The wing is represented as an equivalent plate model, and the Reissner-Mindlin 
displacement field model is used. The Rayleigh-Ritz method is applied to solve the resulting equivalent 
plate problem, with the Legendre polynomials being used as the trial functions. After the stiffness matrix 
and mass matrix are determined by applying the Lagrange’s equations, static analysis can be readily 
performed and the natural frequencies and mode shapes of the wing can be obtained by solving an 
eigenvalue problem. Formulations are such that no limitation is placed on the wing thickness distribution 
as was the case in Cortial 27 . As examples of verifying the present method, a wing-shaped plate, a wing- 
shaped plate with camber, a solid wing and a built-up wing are analyzed respectively and the results are 
compared with those obtained from a detailed FE analysis using MSC/NASTRAN. 
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Assumptions and Formulations 


For the solution of a plate under static or dynamic deformation, the Reissner-Mindlin method, a First- 
order Shear Deformation Theory (FSDT), is based on two assumptions for the displacement field: (1) A 
normal line to the non-deformed middle surface remains a straight line after deformation; and (2) The 
transverse normal stress can be neglected in the constitutive relations. 

According to these assumptions, and assuming linearity, the displacement field of the plate is given as: 


u{x, y, z,t) - u 0 (x , y,t) + z<j> x { x, y,t)' 
v(x, y, Z, t ) = V 0 ( x, y, t ) + z0,(x, y, t ) ► 
w{x,y,z,t) = w 0 (x,y,t) 


(1) 


where as shown in Fig. 1, u,v,w are displacements in the x, y, z direction respectively, subscript 0 refers 
to quantities associated with the plane z = 0 , 0, and 0 y are the rotations about the y and - x axis 


respectively. It is assumed here that the middle surface of the plate is without or with a very small 
curvature, therefore 2 = 0 can be considered to be the middle surface. 

From (1) we can get the strains: 


£ r = — = — 2- + Z- Tx 


du _ du 0 
dx dx 


dv dv 0 

e=-=— Q - + z 


dx 
30 v 


dy dy dy 
dw _ 

r . = 2s . +Z (M. + it) 

" dy dx dy dx dy dx 
dv dw , dw,. 


y =2e = — + — = (/> + 
r>z * dz dy 


3w du 


y dy 

dw n 


( 2 ) 


y — 2e = — h = d> + 

M dx dz dx 


Now we want to analyze a wing by assuming that it behaves as a plate. This assumption is very 
reasonable as long as the wing has a small thickness-chord ratio. For the convenience of calculation, a 
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transformation from (x, >>) to (<*, 77 ) is performed, with the wing configuration in the (x, v) plane, a 
skewed trapezoidal, transformed to a square in the (4,7]) plane, as shown in Fig. 1. 


The transformation can be formulated as 
* = Zw,(£,77)x,. 

y = Xty(£*7)y,- 

1=1 


where 

W 2 (£,/7) = j(l + £)(l-/?) 

N 3 (4,7])=i(l+4)(\ + T]) 

N,(4,7])=i(l-4)(l + 7]) 


Write the inverse of the Jacobian matrix as 


[jy ] 


dx dx 

a? a^ 

dy By 

B4 Br] 


Ju 

J 21 


J 12 
J 22 


( 3 ) 


(4) 


(5) 


We express the displacement components on the plane z = 0 in Eq. ( 1 ), i.e. w 0 , v 0 , w 0 , tf> x and , in 


the following forms 


»o =1 S„ 1 r (<?„) = X£t/,,(<)B,(#)B,(7> 

1=1 7=1 

V„ = ( s„ ) T I = 1 1 v„ (OS, (f )B, 07) 

* = 1 /=! 

“"o = (S„» ) r l = iitV™(0S.(f)S.(/7) • 

m -1 n=l 

!h={B f! ,\ T {q x \=Y^X n w t (Z)Byn) 

p=\ </=l 
J ? 5 

*, = I B„ ) T (?, I = XX (OS, (<f)B, (7) 

r-1 x=l 


(6) 


or 


6 



[u 0 ,v 0 ,w 0 j x ,t y } T =[H][q) 
where I,J,K,L,M,N,P,Q,R and S are integers, 

[q) = [[<lu) T M v ? Mw? Mx\ T MyV? 
{ qu ) = {U u M n ,-,U u ,U 2 l ,-,U 2 J ,-,U n ,-,U,j} T 
[qv) = {Vm-yK L YM»,)={Wu"-W MN ) T 
W x } = {X n ,-,X Pe r,W K } = {F n ,-,F fiJ } r 

is the generalized displacement vector, and 

[m = diag{{B u }\{B KL )\{B MN }\{B PQ } T ,{B ks } t ] 

where 

{B mv } = {B, (£)B, (n), B, (Z)B 2 (tj), ■ • • , B m (Z)B v (77)} t 1 
/iv = U,KL,MN,PQ,RS J 


( 7 ) 


( 8 ) 


(9) 


( 10 ) 


is the Ritz base function vector, in which B,(jc) can be chosen to be the Legendre polynomials: 


B, (*) = />_,(*) 


(ID 


where 



2n + l 

n + l 


P 0 (x) = 1 

P t (x) = x 

xP n {x) —P n _ x (x\n = 1,... 

n + 1 


Strain Energy of a Wing Structure 

The strain energy of a wing structure is 

U =\\\\{G) T {£)dV (12) 

By introducing the stress-strain relations of the generalized Hooke’s law, Eq. (12) becomes 
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u =\\\\{£Y[D]{E}dV 


(13) 


Note that {o } = [D]{£} and [D] r = [D] is assumed, and the integration domain V in Eqs. (12) and (13) 
includes all and only the spaces the components of the wing occupy. 

Using Eqs. (5), we can write 


p/1 


fa/1 


fa/1 




fa/1 
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>=([7] r )-'< 
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a/ 

>=([Jr') T < 

a/ 

► = 
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Ji 
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a/ 
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(14) 


From Eqs. (2) and ( 14), we have 


{£) = 


' ^ ' 


"/n 

J 2! 

0 

0 

0 

0 

zJ u 

zJ 2 \ 

0 

0 

0 

o " 



0 

0 

J 12 

J 22 

0 

0 

0 

0 

zJ \2 

ZJ 22 

0 

0 


- = 

7 12 

J 22 

7„ 
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{e) = [T][e} (15) 


in which {e } is 


f3u 0 3u 0 3 v 0 3v 0 3w 0 3w 0 d<t> x d<p x d<p y ] 

{£) = -T7.T ’TT’T 'Yx'Yy f =l c JWl 

3 £ 3 TJ 3 £ 3 T] 3 £ 3 T} 3 £ dr] d£ drj y \ 


(16) 


where 


[C] = 



0 

0 

0 

0 

0 

[*w 

0 

0 

0 

0 
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0 

0 

0 

0 

0 

[BpQ 'r, ] 

0 

0 

0 

0 

0 


0 

0 

0 
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0 

0 

0 

0 

0 
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(17) 


in which 
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, juv = IJ,KL,MN,PQ,RS, 


l B fjv4rj 1 “ 


'*,'(£)*, 07) 


b;^)b 2 (tj) 

B^B'M) 


B' m (Z)B v (tj) 

B m (Z)B' v (t}) 


{B^ }can be found in Eq. (10), and \q) is the general displacement vector defined in Eq. (8). 

Substitute Eqs. (15), (16) into Eq. (13), and we have 

u = iljj{q) T [CY[TY[D][T][C]{q}dV (18) 

v 

If we write 

U=\{qY[K]{q) (19) 

Then comparison of (18) and ( 1 9) gives 

[K] = JJj [CY [7T [D][T][C]dV (20) 

V 

This is the stiffness matrix of the wing in terms of {q} . Note that the constitutive matrix [D] for different 
parts of the wing structure should be different. 


Kinetic Energy of a Wing Structure 

The kinetic energy of a wing structure is 

T = i JJJ P?dV =i||| pwdw 

v v 


where (v } , the velocity vector, can be written as 


{v} = 



d«0 . Mx 

a t dt 


o 

o 

o 

[**1 

ay 0 a*. 

k 

o 

o 

o 

ri 


dt dt 

dw 0 


_0 0 1 0 0J 


dt 





= [Z][H]{q] 


in which 


( 21 ) 


( 22 ) 
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[Z] = 


(23) 


1 0 0 z 0 

0 1 0 0 z , 

0 0 10 0 


[//] is defined in Eq. (9), and { q } is the time-derivative of { q } 


Then we have 


= ±jjjp{q) T [H] T [ZZ][H]{q}dV 


where 


[ZZ] = [Z] T [Z] = 


1 0 0 z 0 

0 1 0 0 z 

0 0 10 0 


z 0 0 z 2 


0 


0 z 0 0 z l 


Comparing 


T = ±{q} T [M]{q} 

with Eq. (24) we have 
[M] = j\\p[H] T [ZZ][H]dV 


which is the mass matrix of the wing in terms of the general velocity vector, {q} 


(24) 


(25) 


(26) 


(27) 


Numerical Integration of Stiffness and Mass Matrices 

For a specific wing, now we want to evaluate the integrals in Eqs. (20) and (27). 

Assume I-J = K = L- M = N = P = Q- R = S = k, then we know [X - ] and [ M ] will be matrices of 
dimension N x N , where N = 5k 2 . 
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We know coordinates in (x, y) plane are transformed to , but coordinate z remains the same. 


Therefore, for an integral in space (x, y, z ) , we have 
y,z)dV = \' ] \ ] G(Z,?i)d&Ti 

V 


where 


= Xf ! y(Z,rj\z]-\j\dz 


(28) 


(29) 


here N z is the number of integration zones in z -direction, and z^andz^ are integration limits of the i -th 
zone. 

Using the Gaussian quadrature, we can get the numerical value of integral (28) as 


M. N, 






(30) 


,=i H 

where are the weights , ] , rf { ”‘ * are the sampling points, and M g and N g represent the 

number of sampling points used in the t, and 77 directions respectively. 

For a wing composed of skins, spars and ribs, the integrals in Eq. (28) can be detailed as follows: 

(a) Skins 

The sketch of skins at a wing section is shown in Fig. 2. Particulars of integration for skins: 


JJJ F( x, y, z)dV = 1 1 f ■ |s|* + ££ F\j\dz pm 


(31) 


where t LU = t oy [\ + tan 2 a LU , subscript L,U indicate the lower and upper skin respectively. It should be 

noted that for skins composed of laminated layers the skin contributions can be obtained by simply adding 
up the effects of all the layers, with the material constitutive matrix of each layer being treated separately. 
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(b) Spars 


Their contribution to the stiffness and mass matrices can be calculated by performing the integrals for 
each spar. Representative dimensions of a spar are shown in Fig. 3. Thus, for a spar cap, we have: 


JJJ F(x, y, z)dxdydz = J' Y zl ’ M* 




> [ 


(32) 


= \\ dT1 \ l +\ ZV ?' Vfcl&i/ C )f + ^ s ( T lXr}ly[(fJc^ + % s (T)),T}\z\\j\dz 

V-l J- ! \ Zi + * ,L j 

where c is the chord-length at T] : c = ^c 0 (l - 7]) + ^-c, (1 + r}), c 0 : the chord-length at wing root, c, : the 
chord-length at wing tip, and £„(?]) is the spar position function. Equation (32) can be easily computed 

numerically using Eq. (30). 

For the spar web: 


JJJ F(x,y, z)dxdydz = j\ c j_, ^ z] ' \j\d$Ar)dz 

V ^ C Zl+ 2 * L+ 1 

= £ djj L i^x: + ^ ^ ^ + ^ ^ ^ \ j \ dz 


(33) 


(c) Ribs 

The contributions of the ribs to the stiffness and mass matrices can be calculated in a manner similar to 
the one used for spars. The dimensions of a rib are also given in Fig. 3. 

For a rib cap: 


HI F(x, y, z)dxdydz = £ + C-,. V Wf zl W&1& 

V h s\ ' 1 2 / 

= f ^J'Oj/jWJ^T^+J^T 4 * Vfr&»C f 2/ J )7 + J 7r(^)].y&.^2A)7 + 7 7r(^)]*z}-M^ 

J-l J - 1 y 


(34) 


where 5 is the wing semi-span, and T] r (4) is the rib position function. 


For the rib web: 
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( 35 ) 


Jp<* ,y,z)dxdydz = 

= A) 7 + ?,(£)]•?£. (*2 A >7 + *7,(£)1 


The same as for the spars, integration on ribs can be obtained by summing up contributions from all the 
ribs. 


Boundary Conditions and Convergence 

Applying linear and rotational springs with very large magnitudes of stiffness on the boundaries can 
approximate the boundary conditions of clamped edge. Applying linear springs with very large magnitudes 
of stiffness on the boundaries can approximate the boundary conditions of simply supported edge. Details 
of these practices can be found in Lovejoy and Kapania" . 

A series of convergence study on the number of terms, k , of the polynomials used for various cases of 
plates were performed in Lovejoy and Kapania" . The conclusion is, when more terms are used the 
tendency of convergence is obvious, and when k = 8 , the first several modes are very close to 
convergence. In the present research, it is found that there is little difference in the results if the Legender 
or Chebyshev polynomials are used, and the convergence pattern for a built-up wing structure is quite 
similar to that for plates. Therefore it = 8 is used throughout the present study. 


Formulation for Vibration Problem of Wing 

Under the assumption that the wing is a conservative system, the Lagrange’s equations 28 for free 
vibration are 


d_ 

dt 


( dL A 


dL 


dq, I dq, 


= 0 j= 1, - ,//v , fiv =IJ,KL,MN,PQ,RS. 


(36) 
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Since the Lagrangian for the wing is 

L = T - V = T -U (37) 

where V is the potential energy, and by using Eqs. (19), (26) and (36), we can find the natural frequencies 
and mode shapes for the free vibrating wing by solving the following eigenvalue problem 

[K-AM 1jc} = 0 (38) 

where A = a 2 is an eigenvalue of the system of equations, a> is the corresponding frequency in 
radians/second, and { jc } is the corresponding eigenvector. 

Static Problem Solutions 

Assume that an external, distributed force with components P x (x,y,t), P y (x, y,t) and P.(x, 'y,t) is 
applied on the wing structure, then the virtual work done by this load on the infinitesimal area dx dy is 
6Q = 6Q x +6Q y +tQ z (39) 

where 

SQ X = P x (x, y,t) ■ 6u(x, y, z,i) ■ dx ■ dy 

8Q y = P y (x, y, t ) ■ Sv(x, y,z,t) dx ■ dy (40) 

SQ Z =P z (x,y,t)-&v(x,y,z,t)-dxdy 

and by using Eq. (1) we obtain 

<5 Q x = P x ■ (<5 u 0 + z &p x ) dx dy 

5Q y - P y ■ (^ 0 + z ■ &p y ) dx dy (41) 

SQ, = P, Av 0 dx dy 
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Using Eq. (6), that is, approximating the displacements u 0 ,v 0 ,w 0 ,<p x ,<p v in terms of the Ritz functions. 


the total work done by the external force on the whole wing surfaces is given by 




P, {B u ) T [Sq v 1 + P,[B kl ) t {&/ v | + P,{B m ) T 1&, W | 
+ zP,[B rQ ) T {54,} + zP,\B K ) T [fy,] 


nxdy 


= {P) T {Sq) 


where {P} is the generalized load vector 

{/’} = {{/>,} {P 2 ) {*,} W {^ 5 }} r 

in which 

{^*1.2,3 } = JJ ^x.y.z (*’ IJ .KL.MN ) dxdy 
{ n,5 } = ^ JJ P x,y ( X < y> 0 { B PQ.RS 1 7 


(42) 


(43) 


(44) 


If the external force is a concentrated force, the above derivations can be simplified. For instance for the 
first component of the generalized load vector we have 

P, =P x M4 a ^aly^a^a)d){B,j^ a ,ri a )} T (45) 

where (£ a ,ri a ) is the transformed coordinates of the point where the load is applied. 


Using the principle of virtual work, we have, for the static case, the following relation for the 
generalized displacement vector [q] and the generalized load vector { P) 


[*]{<?} = {/>} 


(46) 


Results and Discussion 

In order to assess the accuracy of the present method and test its performance in various situations, a 
series of calculations were carried out for several wing-shaped structures clamped at the root. Results using 
MSC/NASTRAN, a commercial finite element code, are provided for comparison. Finally, for the 
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comparison between the present method and an existing FSDT method, a swept-back box wing used as a 
test case in Livne 21 were calculated for its free vibration and static response analyses. 

Free Vibration Analysis 

(a) A Trapezoidal Plate. The geometric and material parameters for the plate are given as: Span= 192 
inches. Root width=72 inches, Tip width=36 inches. Sweep angle (leading edge)=30°, Thickness=1.8 inch 
(thickness ration at tip= 1/20); Mass Density p = 2.526 x 10~* lb- sec 2 //n 4 . Young’s Modulus 

E = 1 .025 x 1 0 1 lb I in 2 , Poisson's Ratio v = 0.3 . The plate is clamped at the root. 

Comparisons are made in Fig. 4 between the mode shapes as obtained by the present method and those 
by the FEA calculations using MSC/NASTRAN employing 200 shell elements (CQUAD4). The 
comparison of the natural frequencies is also shown in Fig. 4. It can be seen that both the mode shapes and 
natural frequencies as obtained using the present method, are in good agreement with those obtained using 
the FEA. The relative differences of the natural frequencies for the first 8 modes are within -0.62-2. 12%. 

(b) A Trapezoidal Shell with a Camber. All parameters are the same as with the previous case except 
that there is a camber with the camber-chord ratio varying from 2.345% at the root to 0.938% at the tip. 

Comparisons are made between the results as obtained by the present method with those obtained using 
the FEA in MSC/NASTRAN employing 200 shell elements(CQUAD4), as shown in Fig. 5. 

It can be seen that, although this case is very similar to the previous one except for a small camber, 
there are significant differences in the natural frequencies of a number of modes. Most of the variations 
were predicted quite accurately by the present method, as shown clearly in the comparison with the FEA 
results in Fig. 5. But the relative differences were slightly higher than the ones in the previous case, 
varying in a range of - 1 .3 1 % to 5.26% for the first 8 modes. Larger differences for the present case can be 
attributed to the fact that the present method ignores the coupling between the inplane and transverse 
displacements caused by the mid-surface curvature. 
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(c) A Solid Wing. The middle surface of this wing is the same as that of the previous case. Its 
thickness-chord ratio is varied from 0.15 at the root to 0.06 at the tip. The sections were generated by the 
Karman-Trefftz transformation 29 . 

Comparisons are made in Fig. 6 between the results as obtained by the present method with those 
obtained by the FEA calculations using MSC/NASTRAN employing 250 solid elements (CHEXA and 
CPENTA) and 572 nodes. It can be seen that although there are thickness variations as well as a camber, 
the present method yields results that compare quite well with those obtained using the FEA. The relative 
differences for the first 8 modes were within -5.82-1.42%, comparable to those in the previous case. 

(d) A Built-up Wing Composed of Skins, Spars and Ribs. The outside geometrical shape is the same as 
in the previous case, the solid wing. There are 4 spars and 10 ribs distributed uniformly in the wing. 
Particulars of the wing are: Skin Thickness / 0 = . 1 18m ; Spar Cap Height /i, = . 197m , Spar Cap Width 

/, = .373 in , Spar Web Thickness t, = .058m ; the ribs have the same cap dimensions and web thickness as 
the spars. 

The FEA calculations are made by using MSC/NASTRAN employing 370 elements and 1 10 nodes. 
The wing skins were modeled using shell elements (CQUAD4), the spar and rib caps were modeled using 
bar elements (CBAR), and the spar and rib webs were modeled using shear panel elements (CSHEAR). 
Comparison between the mode shapes as well as the corresponding natural frequencies as obtained by the 
two methods are shown in Fig. 7. It can be seen that the mode shapes were simulated equally well by the 
present method as compared to the FEA, and it is found that the relative differences for the first 8 modes 
were within -4.79-2. 15 %, comparable to those in the previous cases. 

(e) A Box Wing used as a test case in Livne 21 . This is a cantilevered ail aluminum wing swept back by 
30° . It has a constant thickness and constant chord length. Its 5 spars and 3 ribs with identical cross 
sections are bonded to the top and bottom cover skins. Details of this box wing can be found in Refs. 2 1 
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and 30. The same kinds of elements were employed as in the previous case. Results for the natural 
frequencies by the FEA using MSC/NASTRAN and the present FSDT are shown in Table 1 , in 
comparison with those given in Livne 21 by the FEA using ELFINI and a FSDT based on simple- 
polynomial trial functions. While there are some differences between the two FEA calculations, which may 
have been caused by different discretization and element choices, the accuracy of the present FSDT results 
are promising. 

Displacement under Static Loads 

The built-up wing in case (c) of the free vibration analysis is used here. Three cases of static load were 
considered and the results compared with the FEA calculations using MSC/NASTRAN are shown as 
follows. 

(a) Tip Point Force. A downward (- z -direction) force of magnitude of 1 lb is applied at the middle 
point of the wing tip. The displacements along the leading and trailing edge of the wing are shown in Fig. 

8. It can be seen that the present method calculated the vertical displacement w accurately compared with 
the FEA, and also predicted quite well the trends of variation for the other two displacement components, 
u and v. 

(b) A Force Distribution. A downward (- z -direction) force of magnitude of 1 lb is applied at every 
upper-surface nodes of the FEA model. This is a case similar to the wing being under uniform pressure 
difference between its upper and lower surfaces. The displacements along the leading and trailing edge of 
the wing are shown in Fig. 9. Quite similar results to those in load case (a) were obtained. 

(c) Tip Torque. A downward (- z -direction) force of magnitude of 1 lb is applied at the tip of the fore- 
most spar, while an upward ( z -direction) force of magnitude of 1 lb is applied at the tip of the aft-most 
spar. This is a case in which the wing tip is subjected to a torque. The displacements along the leading and 
trailing edges of the wing are shown in Fig. 10. The relative difference for w at the tip/leading-edge comer 
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and tip/trailing-edge comer are 26.4% and 2.64% respectively. Note that the absolute differences are 
0.397 x 10" 4 inch and 0. 1 8 1 x 10 -4 inch respectively; therefore the large relative difference at the 
tip/leading-edge comer is because of the small magnitude of w along the leading edge, and the difference 
between the twist angles predicted by the two methods would be small. 

(d) The Box Wing in Livne 21 . The swept-back box wing is under a downward point force of 1 lb at the 
tip of the rear spar. Displacements at the wing tip from measurement, FEA calculation and present method 
are shown in Table 2. It can be seen that the present method yielded very good results for this test case. 
Skin Stress Distributions 

The upper and lower skin stress of the wing in (a) of the above static cases were calculated using the 
present method. The Von Mises stress distribution along a line with a distance of 5% span to the root chord 
is shown in Fig. 1 1(a) in comparison with points obtained using MSC/NASTRAN. Also the Von Mises 
stress distribution along a span-wise line with a distance of 37.5% chord length to the leading edge is 
shown in Fig. 1 1(b) in comparison with points obtained using MSC/NASTRAN. It can be seen that, 
although there are substantial differences (the largest one is about 15%) between the present calculations 
and the stresses determined using the FEA, the variation trends of the stresses from both methods are quite 
similar. This means that the position of the largest stress determined by the present method will be reliable. 

Conclusion 

A method capable of static and vibration analyses of the general built-up wing stmctures composed of 
skins, spars and ribs has been developed and comparisons for a series of examples with commercial FEA 
calculations have shown the accuracy of the method for design purposes. On the assumption that the wing 
structure behaves like a plate whose deformation can be modeled by the FSDT of Reissner-Mindlin, the 
Rayleigh-Ritz method is applied to solve the plate problem, with the Legendre polynomials being used as 
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the trial functions. The stiffness and mass matrices are determined by applying the Lagrange’s equations, 
and can be calculated numerically by using the Gaussian integration quadrature. Then static analysis can be 
readily performed and the natural frequencies and the mode shapes of the wing can be obtained by solving 
an eigenvalue problem. 

Comparison of results by the present method with those by the commercial finite element analysis code 
MSC/NASTRAN for a series of 5 vibration problems, 4 static loading problems, and 1 stress distribution 
problem showed an overall good agreement between the two approaches with different methodologies. 
Mode shapes and natural frequencies for cases from a thick wing-shaped plate, the same plate with a 
camber, a solid wing, to built-up wing structures composed of skins, spars and ribs, have all shown that the 
present method has a fairly good correlation to the FEA, although results for simpler cases seem to be 
more accurate. It was also shown that static displacements and stress variation trends of wing structure can 
be predicted by the present method quite accurately. The present method is formulated mostly in matrix 
form and calculation can be readily carried out in the MATLAB environment. It is suitable to be used for 
the early stages of wing design. Due to the efficiency of the method, it can also be used as a means to 
analyze the shape sensitivity of wing structures. 
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Table 1 Natural frequencies (Hz) of the cantilevered swept-back box wing 


Mode 

Description of 

FEA 

FSDT 

FEA 

FSDT 

No. 

Mode Shape 

(Livne 21 ) 

(Livne 21 ) 

(present) 

(present) 

1 

I s ' bending 

115.6 

114.7 

116.6 

118.0 

2 

In plane 

317.6 

312.4 

327.9 

349.7 

3 

l sl torsion 

418.4 

428.9 

409.4 

419.1 

4 

2 nd bending 

576.4 

575.3 

572.1 

571.1 

5 

2 nd torsion 

1086 

1125 

1064" 

1090 


Table 2 Displacement (in) of the cantilevered swept-back box wing 


Position 

Measured 21,30 

FEA(present) 

FSDT(present) 

Front spar tip 

1.81 x 10"* 

1.79x10^ 

1.74x10“* 

Rear spar tip 

2.21x10^ 

2.19x10“' 

2.20x10“' 




Fig. 1 The Coordinate System and Its Transformation 
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Mode No. 1 /=! 0.449 Mode No. 1 /=10.461 



by FEA 


Equivalent Plate Model 


Fig. 4 Mode Shapes and Natural Frequency / {rad! s) for a Trapezoidal Plate 






Mode No. 3 / *97.068 





by FEA 


Equivalent Plate Model 


Fig. 5 Mode Shapes and Natural Frequency / (rad / s) 
for Wing-Shaped Shell with a Camber 
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by FEA Equivalent Plate Model 


Fig. 6 Mode Shapes and Natural Frequency / ( rad / s) for the Solid Wing 
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Mode No. 1 /= 54.2 














by FEA 


Equivalent Plate Model 


Fig. 7 Mode Shapes and Natural Frequency / ( rad / s) 
for a Built-up Wing Composed of Skins, Spars and Ribs 
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1 0Ou, 1 0Ov, w (inch) 1 00u - 1 00v - w ( lnch ) 


-3 

x io 



-3 

xIO 



Fig. 8 Comparison of Displacements for Load Case of Tip Point Force 
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0.01 



(a) leading edge 



Fig. 9 Comparison of Displacements for Load Case of a Force Distribution 


3 





(a) leading edge 



Y (inch) 

(b) trailing edge 


Fig. 10 Comparison of Displacements for Load Case of Tip Torque 
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(a) Near the Root Chord 



(b) Near the Central Spar 


Fig. 11 Comparison of Von Mises Stress on the Upper and Lower Skins 
of a Wing under a Point Force at the Wing Tip 
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Fig. 1 The Coordinate System and Its Transformation 

Fig. 2 Wing Skin 

Fig. 3 Wing Spar or Rib 

Fig. 4 Mode Shapes and Natural Frequency / 

Fig. 5 Mode Shapes and Natural Frequency / 

Fig. 6 Mode Shapes and Natural Frequency / 

Fig. 7 Mode Shapes and Natural Frequency / 

Spars and Ribs 

Fig. 8 Comparison of Displacements for Load Case of Tip Point Force 
Fig. 9 Comparison of Displacements for Load Case of a Force Distribution 
Fig. 10 Comparison of Displacements for Load Case of Tip Torque 

Fig. 11 Comparison of Von Mises Stress on the Upper and Lower Skins of a Wing under a Point 
Force at the Wing Tip 


{rad / s) for a Trapezoidal Plate 

( rad / s ) for Wing-Shaped Shell with a Camber 

{rad / 5 ) for the Solid Wing 

{rad Is) for a Built-up Wing Composed of Skin, 
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